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(54) Determining height information 

(57) A method is provided for phase unwrapping in 
coherent imaging systems which use two complex 
images to form an interferogram. In the method, the 
interferogram of phase differences is divided into a rela- 
tively large plurality of transform blocks, wherein each 
transform block comprises a matrix of o x p data sam- 
ples, with o and p = 2 n + 1, and n being an integer 
greater than 2. Further, adjacent transform blocks are 
preferably defined to partially overlap. The wrapped 
phase values of adjacent data samples within each 
transform block are then compared and if the difference 
therebetween exceeds a predetermined value, the 
entire block of data samples may be discarded for fur- 
ther use. Data samples of the retained transform blocks 
are then unwrapped to obtain phase values via an 
unweighted, least-squares technique implemented by 
Fast Fourier Transform. Using a known height value cor- 
responding with a single data sample, an integration 
constant can be determined for a corresponding first 
transform block. Path-following may then be used to 
determine an integration constant for the other retained 
transform blocks. The phase values and integration con- 
stants may then be employed to determine height val- 
ues for the interferogram. The method yields enhanced 
accuracy with reduced computational burden. 
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Description 

[0001] The present invention relates to a method and system for determining height information. An embodiment of 
the present invention relates to coherent imaging systems in which two complex images of an imaged region may be 

5 merged to interfere in such a way as to cancel the image components which are common to both images and recover 
the information that is transduced by image-domain phase data. More particularly, the embodiment pertains to an 
improved method for phase unwrapping in such imaging systems, and is particularly apt for application to interferomet- 
ric synthetic aperture radar (IFSAR) applications to determine terrain height within an imaged terrain region. In this 
regard, the embodiment implemented to supplement and thereby improve known IFSAR processing systems, such as 

10 the systems disclosed in C. Jakowatz, Jr., D. Wahl. P. Eichel, D. Ghiglia and R Thompson, SPOTLIGHT -- MODE SYN- 
THETIC APERTURE RADAR: A Signal Processing Approach (1996), and D. GHIGLIA and M. PRITT. TWO DIMEN- 
SIONAL PHASE UNWRAPPING THEORY, ALGORITHMS AND SOFTWARE (1998). 

[0002] Computed imaging systems are now widely employed in a variety of applications, including medicine, astron- 
omy and terrain analysis. Imaging modalities used for such applications include computer tomography, magnetic reso- 
15 nance imaging, ultrasound, synthetic aperture radar and radio astronomy. 

[0003] Of particular interest here, synthetic aperture radars (SARs) have been employed to produce high-quality 
images of the earth's terrain. Such images are obtained by overhead transmission/ receipt of pulses of microwave 
energy at a predetermined frequency. In this regard SARs provide particular functionality due to their capability to 
image in darkness and to penetrate cloud-cover. Further, because SAR systems use a form of coherent illumination. 
20 SAR systems are capable of transducing the complex reflectivity of terrain within an imaged region. In such applica- 
tions, the reflectivity function is modulated by phase terms that are dependent upon the imaging system geometry. As 
a result, when two SAR images are made of an imaged region, it is possible to interfere the two complex images in such 
a way as to cancel the scene reflectivity which is common to both images and recover the information that contains the 
scene topography transduced by the image-domain phase data. Such systems may be generally referred to as interfer- 
es ometric synthetic aperture radar (IFSAR) systems. 

[0004] IFSAR systems, both aircraft and space borne, have been used with moderate success to date to provide ter- 
rain height for regions on the earth's surface. Such systems may consist of a single vehicle having one radar transmitter 
and two spaced receive antennas mounted thereupon, wherein two complex images of an image region may be 
obtained upon a single pass of the vehicle over an imaged region. Alternatively, the system may comprise one or more 
30 .vehicles that each have one radar transmitter and one receive antenna mounted thereupon, wherein two complex 
images are obtained by passing over an imaged region twice. In either case, two complex images of the same region 
are formed. After acquisition, the images may be registered such that the phase differences between corresponding 
image pixels, or data samples, may be extracted to form an interferogram. As will be appreciated, the phase differences 
reflected by the interferogram are wrapped. That is, the phase differences are ambiguous module two pi (2n). 
35 [0005] In order to derive height information from the interferogram, the wrapped phase differences must be 
unwrapped and corresponding integration constants must be determined. As such, phase unwrapping should be com- 
pleted in a manner that resolves the 2n ambiguities so that unambiguous terrain heights can be assigned to the phase 
values. In addressing such task, it has been recognized that the imaged terrain cannot be of a nature that yields phase 
values that exceed the Nyquist rate without adversely impacting the accuracy of results. Specifically, adjacent sample- 
40 to-sample phase differences of unwrapped interferometric data should be no more than 180°. Such limitation can 
become problematic when the imaged region comprises steep pastoral terrain (e.g., near vertical natural geographical 
features) or cultural features (e.g., man-made structures such as buildings). When such features are present, phase 
unwrapping may result in inconsistent data that renders the entire height estimation unreliable. 

[0006] To understand such inconsistencies, consider a closed path consisting of one step forward, a step to the left, 
45 a step to the left again, and then a f inal step to the left. After the four steps, one should arrive back at the starting point. 
Interferometric differential phase data is supposed to represent terrain height, but in situations that present the above- 
mentioned problem (i.e., adjacent samples whose phase difference is >180 °), it is possible that the sum of phase dif- 
ferences around a closed, four-point path in the interferogram is non-zero. Such a result would indicate that if one con- 
verted the phase differences to height differences and summed them around the path, one would not arrive back at the 
so starting height. 

[0007] Existing phase unwrapping algorithms are of two general types: least squares and path following. Least- 
squares algorithms determine the phase surface which best fits the ensemble of pixel-to-pixel phase differences over 
the entire interferogram. If inconsistencies of the above-noted nature are present, the least-squares process attempts 
to minimize their deleterious effects by minimizing the residual fitting error. Path-following algorithms, on the other hand, 
55 numerically integrate the pixel-to-pixel phase differences over the interferogram. in the process either avoiding or mini- 
mizing inconsistencies by selecting closed paths where error is minimized. 

[0008] Systems based on these methods have not been able to meet the requirements of many potential applications. 
To date. IFSAR systems have achieved accuracies of a few meters on the average. But because of "errors" in the 
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unwrapping algorithms (i.e., due to the noted inconsistencies), accuracies may be excellent in one region of the image 
and quite poor in another, and one has no way of knowing which regions are good and which are bad. Further, future 
systems are desired which can provide sub-meter accuracies with a high degree of assurance. The noted phase 
unwrapping techniques are not up to the task 

5 [0009] In this regard, both least-squares and path following algorithms can be undermined by terrain characteristics 
that are relatively common. By way of example, if the number of inconsistencies in an imaged region is large, the least- 
squares solution is impractically crude because the effects of the bad data are smeared throughout the image. If one 
has a priori knowledge of regions of defective data, then a weighted least-squares can be used which assigns low, or 
zero weight to poor data. To date, however, no one has developed a robust way to determine appropriate weights from 

10 the interferogram. 

[0010] Path-following methods are also undermined by large numbers of inconsistencies because automated tech- 
niques for finding satisfactory integration paths fail, and the underlying algorithm cannot complete the unwrapping proc- 
ess. One does not know in advance when such failures will occur. 

[0011] A final difficulty with existing phase unwrapping methods relates to the need to incorporate large numbers of 
15 tie points (i.e., locations within the imaged region having known heights) into the algorithm employed in the correspond- 
ing systems. To ensure the best possible accuracy, one should incorporate as many known height tie points as possible. 
Neither least-squares nor path-following systems lend themselves to incorporating dense grids of tie points, and this will 
be needed in future systems which aim for sub-meter accuracies. 

[0012] Aspects of the invention are specified in the claims. Another aspect of the invention provides a system 
20 arranged to carry out the method of claim 1 and its dependent claims. 

[001 3] An embodiment of the present invention seeks to provide an improved method for phase unwrapping in imag- 
ing systems that employ two complex images, including in particular interferometric synthetic aperture radar (IFSAR) 
systems. In particular, a primary objective of the present invention is to provide for improved accuracy in height deter- 
minations in IFSAR systems. 

25 [0014] A further embodiment of the present invention seeks to provide a method of phase unwrapping that reduces 
the use of inconsistent data, thereby yielding enhanced accuracy (e.g., in height determinations for IFSAR applica- 
tions), 

[0015] An additional embodiment of the present invention seeks to provide a phase unwrapping that can be carried 
out in a computationally efficient manner. 

30 [001 6] Yet a further embodiment of the present invention seeks to provide a method of phase unwrapping that allows 
the use of multiple tie points in a convenient and effective manner. 
[001 7] An embodiment of the present invention comprises a method which 

facilitates the extraction of dimensional, or positional, information from an interferogram generated from two 
complex images of an imaged region. The method includes the step of dividing the interferogram into a plurality of sub- 

35 regions, or transform blocks, wherein each transform block comprises a plurality of complex data samples that each 
include a wrapped phase value. In this regard, it has been determined that each transform block should comprise an o 
x p matrix of data samples, wherein o and p should be equal to 2 n + 1 , wherein n is an integer preferably greater than 
2, and most preferably greater than 3. Further, each transform block may be advantageously defined to partially overlap 
each of its neighboring, or adjacent, transform blocks. In this regard, the overlap regions should preferably comprise a 

40 number of data supplies that is between about 5% and 50% of the total number of data samples comprising each trans- 
form block. As will become apparent, the use of a transform block approach to phase-unwrapping yields a number of 
benefits in the present invention. Such benefits are particularly apt for use in IFSAR systems, and the present invention 
will be further summarized in relation thereto. 

[0018] In one aspect of the present invention, the method comprises unwrapping a phase difference value i) for each 
45 of the plurality of complex IFSAR data samples comprising at least a first transform block of an IFSAR interferogram. 
including an overlapping region with at least a second transform block of an IFSAR interferogram, and it) for each of the 
plurality of complex data samples comprising at least the second transform block, including the noted overlapping 
region. A known height value, or tie point, corresponding with one of the plurality of complex data samples comprising 
the first transform block is then used to establish an integration constant for the first transform block. Such integration 
so constant can then be employed with the unwrapped phase difference values for the overlapping region to determine an 
integration constant for the second transform block (i.e., since the average of the unwrapped phase values for the over- 
lapping region ot the first transform block can be assumed to be equal to the average of the unwrapped phase values 
for the corresponding overlapping region of the second transform block). Using this methodology, an integration con- 
stant for other ones of the plurality of transform blocks comprising an interferogram can be determined via path-follow- 
55 ing. As will be appreciated, a height value corresponding with a given complex data sample can be determined using a 
corresponding unwrapped phase value and an integration constant for the transform block within which the sample is 
located. 

[001 9] In another aspect of the present invention, the method comprises the step of analyzing each of the plurality of 
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transform blocks comprising an IFSAR interferogram to identify inconsistencies, or ambiguities, in the corresponding 
data samples. More particularly, such analyzing step may comprise the following substeps that may be carried out prior 
to any phase unwrapping of the complex data samples: i) determining a difference between the wrapped phase values 
of each set of adjacent complex data samples comprising a given transform block, and //) comparing each of the differ- 

5 ence values to a predetermined value, wherein all complex data samples comprising a transform block are discarded 
from further use in the method when any single difference value for the transform block exceeds the predetermined 
value, and wherein all complex data samples comprising a transform block are retained for further use in the method 
when all difference values corresponding with the transform block are less than the predetermined value. Such prede- 
termined value may be preferably set to be less than about 1 35°. and even more preferably less than about 90°. As will 

to be appreciated, by discarding transform blocks comprising inconsistent data in the described method, the accuracy 
achievable by the present invention is enhanced. 

[0020] In a related aspect of the present invention, the method may comprise the step of unwrapping phase difference 
values for the data samples comprising each retained transform block (e.g. data samples that have not been discarded 
pursuant to the above-noted analyzing step), via use of an unweighted least-squares algorithm. Preferably, such algo- 
15 rithm may be implemented by Fast Fourier Transform. As a related consideration, the total number of data samples 
comprising each transform block should preferably comprise a robust statistical sample of interferometric phase differ- 
ences (e.g., preferably more than about 100 samples, and most preferably more than about 200 samples). On the other 
hand, transform blocks should not be excessively large, so that when ambiguities are found, the discarded region is not 
excessively large. 

20 [0021 ] In one application, an IFSAR system is employed to acquire two complex images of an imaged terrain region. 
By way of example, such system may comprise a single transmitter for transmitting microwave energy pulses from an 
airborne or space-borne vehicle toward an imaged region, and a single receiver for receiving microwave energy 
reflected from within the image terrain region. In such an arrangement, two passes over the imaged region would be 
necessary to obtain the two complex images. Alternatively, two complex images may be obtained via a single pass 

25 when a single transmitter with two receivers are used on the vehicle. After acquisition of the two complex images, the 
images are spatially registered so that complex data samples corresponding with a given location in the imaged terrain 
region can be further used. More particularly, following registration, an interferogram can be generated by merging the 
two complex images. As will be appreciated, the noted complex image merging serves to extract wrapped phase differ- 
ences between the corresponding pixels or data samples comprising the two complex images. 

30 [0022] To facilitate processing, the interferogram may be smoothed using a scanning filter. The interferogram is then 
divided into a plurality of partially overlapping transfer blocks. By way of example, where an interferogram comprises a 
matrix of data samples, or pixels, of at least about 4,000 by 4,000, the transform blocks may each be 1 7 x 17, or 33 x 
33 data samples. Additionally, the overlapping regions may be on the order of 4 data samples wide. 
[0023] Upon division of the IFSAR interferogram into the overlapping transform blocks, each transform block may be 

35 analyzed for purposes of identifying data inconsistencies (e.g. arising due to data sampling from terrain that exceeds 
the Nyquist rate). In this regard, the wrapped phase difference values comprising each set of adjacent data samples of 
the transform block may be compared to determine if the difference therebetween exceeds a predetermined value. In 
the event any given difference value exceeds the predetermined value (i.e. thereby indicating inconsistent data), the 
entire transform block of data samples is discarded from further use in the method. In this regard, it has been found that 

40 a predetermined value of about 135° or less, or even about 90° or less may be advantageously employed in such anal- 
ysis. Upon discarding the transform blocks comprising inconsistent data, phase unwrapping of the data samples com- 
prising the retained transform blocks may be conducted using an unweighted least-squares algorithm, as implemented 
by Fast Fourier Transform. 

[0024] After unwrapping, phase integration constants corresponding with each retained transform block may be deter- 
45 mined. More particularly, in the described IFSAR application, such integration constants may be determined via use of 
a single-known tie point (i.e., a known height corresponding with a single complex data sample) in combination with a 
path-tollowing approach. In this regard, it is again noted that since each transform block overlaps its neighbor, the 
avenge unwrapped phase values for the corresponding overlap region of any two adjacent transform blocks should be 
equal. Using this fact, together with the integration constant established using the single known tie-point, the integration 
so constants for other transform blocks may be determined. 

[0025] In the event that a plurality of adjacent transform blocks of data samples have been discarded (i.e., in the ana- 
lyzing step) so as to form a band between different regions, or sets, of retained transform blocks, it can be appreciated 
that one or more sets of retained transform blocks may be rendered "inaccessible" by path-following. In such a situation, 
a further single known tie point located within an "inaccessible" set of transform blocks may be used to determine the 
55 integration constant for the corresponding transform block, and in turn, integration constants may be determined via 
path-following for the other transform blocks comprising the "inaccessible" region. If a tie point is not known for a given 
"inaccessible" region, the transform blocks and related data samples comprising such may be discarded (i.e., not 
retained) from further use in the method. 
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[0026J After determination of the integration constants for each of the retained transform blocks, the unwrapped phase 
values (in radians or degrees) and integration constants may be used together with a pre-determined phase-to-height 
conversion factor to determine a height value for each data sample comprising the retained transform blocks. The final 
product is a standard Digital Elevation Model (DEM) of the terrain. 

5 [0027] It should be noted that multiple tie points may be readily used at this point in the process to further enhance 
accuracy. That is, where more than one tie point is known for a set, all adjoining, returned transform blocks, the known 
tie points that were not previously used (i.e., in the step of determining integration constants) now may be advanta- 
geously employed. The additional tie points may be incorporated into the DEM by constrained least-squares methods 
which adjust the transform block heights (i.e., dimensions) to reduce, or minimize, the squared height differences in the 

10 overlap regions. The least-squares constraints are the known tie points. 

[0028] Numerous extensions, additions and advantages of the present invention will become apparent to those skilled 
upon further consideration of the illustrative description that follows with reference to the accompanying drawings, in 
which: 

is Fig. 1 is a process flow diagram of one embodiment of the present invention. 

Figs. 2A and 2B illustrate an exemplary interferogram and a portion thereof divided into overlapping transform 
blocks, respectively. 

Fig. 3 illustrates process flow substeps of an alternate approach for step 80 of the embodiment of Fig. 1 . 
Fig. 4 illustrates process flow substeps for step 90 of the embodiment of Fig. 1. 

20 

[0029] Figs. 1 -4 are directed to an interferometric synthetic aperture radar (IFSAR) system embodiment of the present 
invention. As will be appreciated, other applications of the present invention may include magnetic resonance imaging 
(MRI) and astronomical imaging. 

[0030] With reference to Fig. 1 , the illustrated IFSAR process embodiment 10 comprises a number of steps. Initially. 

25 an IFSAR image acquisition system is used to acquire two complex images of an imaged terrain region (step 20). In this 
regard, and as will be appreciated, either aircraft or space-borne vehicles may be used to acquire the two complex 
images. The images are acquired by transmitting microwave energy pulses at a predetermined frequency toward the 
imaged terrain region, and receiving resultant microwave energy reflected from the terrain within the region. For each 
complex image, the received energy is detected in spatial relation to the to the imaged region, demodulated and stored 

30 to yield a two-dimensional array of complex sample data. By way of example, a single vehicle may be used with one 
microwave pulse, or radar, transmitter and two spaced receive antennas to acquire the two complex images during a 
single pass over the imaged region. Alternatively, one or more vehicles may be employed that each have a radar trans- 
mitter and one receive antenna may be employed to obtain two complex images by passing over the imaged region 
twice. After image data acquisition, demodulation and storage, an image-formation processor may be employed to com- 

35 plete the processing steps contemplated by the illustrated embodiment 10. Such processing may be completed either 
onboard an imaging vehicle, at a ground-based location, or at some other location remote from the imaging vehicle. 
[0031 ] As will be appreciated, the described IFSAR image acquisition system uses a form of coherent illumination to 
transduce the complex reflectivity of the imaged terrain region. Such reflectivity function is modulated by phase terms 
that are dependent upon the imaging system geometry. As a result, when two complex images are acquired by the 

40 IFSAR system, the two complex images may be merged to interfere in such a way as to cancel the scene reflectivity 
which is common to both and to recover the geometric information that contains the topography of the imaged region 
as transduced by the image-domain phase data. Such phase data is employed to determine height values within the 
imaged region. 

[0032] In this regard, following acquisition of the two complex images (step 20), the embodiment 10 of Fig. 1 provides 
45 for 2D registration of the images (step 30). This may be done in a variety of ways known in the art. For example, regis- 
tration may entail the generation of a set of control points or local image-to-image displacement vectors, the calculation 
of a warping function, and image resampling. In this regard, registration control points may be generated via two-dimen- 
sional correlation of image subregions. Warping functions employed in the registration step may be polynomial-based, 
spline-based, or a combination of the two. A simple bilinear interpolator may be employed for resampling the source 
so image(s) according to the warping function. 

[0033] In conjunction with registration, step 30 further comprises the formation of complex difference samples, which 
consist of multiplying each complex pixel from one image by the complex conjugate of the associated pixel from the 
other image. Such formation of complex difference samples may be presented as follows: 

55 0 iJ = Arg[P iJ • Q tf ] 

where: 
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Pj j= sample from one image 

Q j j = complex conjugate of sample from other image; and 

where the matrix, or plot of 0 j j determines an interferogram. 

5 [0034] The real and imaginary components of the complex difference samples may then be smoothed (step 40) by a 
moving averager with dimensions 3 by 3, 5 by 5, or 7 by 7 sample, depending on the quality of the data. 
[0035] As shown in Fig. 1, the process embodiment 10 further comprises the important step of dividing the interfero- 
gram into a plurality of partially overlapping transform blocks (step 60). In this regard, Fig. 2A illustrates an exemplary 
interferogram 61 comprised ot a plurality of pixels, or complex difference data samples. By way of example, interfero- 

10 gram 61 may comprise anMxN matrix of data samples, wherein M and N = 4000. For purposes of further explanation, 
a portion 62 of the interferogram 61 has been enlarged in Fig. 2B. 

[0036] Fig. 2B includes a corner region 63 of the interferogram portion 62 with a plurality of overlapping transform 
blocks 64 defined therein. As can be seen, the transform blocks 64 are overlapped to define overlapping regions 65. In 
one arrangement (e.g., wherein M and N = 4000), each transform block 64 may comprise m x n pixels, or complex data 
is samples, (e.g., wherein m and n = 17). In such arrangement each transform block 64 may be established to overlap 
each of its neighbors by a predetermined band of pixel rows and columns, e.g., each predetermined band being 4 pix- 
els, or data samples, wide. 

[0037] As will be appreciated upon consideration of the description that follows, the decomposition of interferogram 
61 into a relatively large number of overlapping transform blocks 64 yields significant processing efficiencies and other 

20 advantages. In this regard, the total number of pixels, or data samples, in each transform block 64 should preferably be 
at least 100 samples. Further, the total number of pixels, or data samples, in overlapping regions 65, should preferably 
represent between about 5% to about 50% of the total number of pixels in a transform block 64. It should be noted that 
while equal-sized transform blocks and equal-sized overlapping regions yield certain processing advantages, the 
present invention may be employed with differing sizes as well. 

25 [0038] Upon defining overlapping transform blocks in an interferogram (e.g., 61 ), the next step in process embodiment 
10 is to analyze each transform block (e.g., 64) to identify "unacceptable" inconsistencies, wherein blocks which contain 
such inconsistencies are discarded from further analysis or use in the process (step 70). To conduct such analysis, all 
adjacent (i.e., sample-to-sample) first-order, wrapped phase differences in an interferogram (e.g., 61) are computed. 
Such first-order wrapped phase differences are ten incfividually assessed to determine whether or not they lie outside 

30 of a predetermined acceptance range. In this regard, it will be appreciated that if each data sample has been provided 
at the Nyquist rate with respect to the spatial frequency content of the terrain, the legitimate range of sample-to-sample 
differences should be between -180° and 180°. In this regard, however, it has been recognized that higher frequency 
differences between adjacent data samples may alias down into the 180° to 180° range and become indistinguishable. 
To reduce the possibility of using such contaminated data, the above-noted, specified acceptable range for use in 

35 embodiment 1 0 may be defined as follows: 



-T£ A 0 £ T°; 



wherein T is preferably about 135° (e.g., for 33-by-33 point transform block), and wherein T is even more preferably 
40 about 90° (e.g., for a 1 7-by-l 7 transform block). If the wrapped first-order phase difference A 0 between any two adja- 
cent samples in a given transform block fall outside of the above-noted acceptable range, the data associated with all 
of the samples in the associated transform block (e.g.. 64) is discarded from further analysis. As will be apparent to 
those skilled in the art, the described approach for handling inconsistencies in the process embodiment 10 yields 
enhanced accuracy in height estimations. 
45 [0039] Referring again to Fig. 1 , the process embodiment 10 comprises the further step of unwrapping the phase for 
each pixel, or data sample, within each retained transform block (step 80) of an interferogram (e.g., 61). In this regard, 
and particularly since inconsistent data has been discarded, phase unwrapping may be advantageously completed 
using an unweighted, least-squares procedure, implemented by Fast Fourier Transform (FFT). As will be appreciated, 
the combination of relatively small transform blocks (e.g., 64) in an FFT implementation results in a relatively low corn- 
so putational burden for completing the unwrapping step 80 of the process embodiment 1 0. 

[0040] More particularly, in completing the phase unwrapping step 80. an algorithm may be employed that is based 
on a Z-transform solution of the difference equation that approximates the phase 0 at each sample data point in terms 
of its 4 neighbors. From analytic geometry the slope of the line between points (x1 ,y1) and (x2,y2) in an interferogram 
(e.g.. 61) may be defined as follows: 



55 
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This lets us predict other points on the line from: 

y-yi= m i.2-( x - x i)- 

[0041] Now, consider an array of five points in a plane of wrapped phase values fw(j,k). Using the point-slope equa- 
tion, we can make four independent predictions of the unwrapped phase value <t>G,k) based on its four nearest neigh- 
bors: 

* M=4>M- 1 * fw jX fw j,k-J=*j,k- 1 -"*jhm ( 3 ) 
Where the first order differences are defined as: 

Under general conditions, the maximum-likelihood estimator of the expected value of a population based on a finite set 
of samples is the ordinary average. Accordingly, if we average the four independent estimates we get the combined esti- 
mator: 

25 4 -*M = */-1,fc + +/ + 1.fc + *y,/t-1 + +/,/t + 1 # ' * + m jj,k- m J/-\k +mk j,k' mk j,k-1 ■ ( 6 > 

[0042] This is precisely the difference equation which must be solved to obtain an unweighted least-squares estimate. 
The least-squares formulation determines the estimates <t>G.k) such that the sum of the squared differences between 
estimated and measured phase slopes is minimized. 
30 [0043] The Z-transform is widely used as a way to solve difference equations. The above-noted difference equation 
is noncausal in that <Kj,k) is given in terms of both left-hand and right-hand values <KM,k) and <Ki+1,k). Hence a two- 
dimensional, two-sided Z-transform is appropriate. The definition is Z = F(z,u) where: 

OO CO 

35 F{z,u)= £ £ fwjy z • u* . (7) 

y=-t» ft =-oo 

[0044] Anticipating a solution using the Discrete Fourier Transform (DET), we transform to one-sided summations with 
40 the result: 
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F<*u)-£[ £ [{fw jfk • u ~ k +fw jrk • u *) • z ' + (fw , itk • a ' k + fw ^ - u * ) • z ']! 

y=oU=0 J 



(8) 



45 



where we have ignored the repeated points at j = 0 and k = 0. To force a real, periodic solution, we require mirror-image 
symmetry in the input domain so that the Z-transform is real on the unit circle. To create this symmetry, the data is rep- 
licated over the other three quadrants, doubling the dimensions of the array. After this expansion of the data, we have 
so fG,k)=f(-j,k)=fQ,-k)=f(-j,-k) , and the above summation becomes: 

OO OO 

F(z,u)= £ £ fw jy [(u -*+u *) • (z + z')] (9) 

y=o/c=o 

55 

[0045] After the mirror-image extension, F(z,u) is real and even when evaluated on the unit circle. This is because 
fQ,k) is real and even, and the z and u terms become real cosines on the unit circle (e.g., z = expfltoAt ). Continuing with 



7 



ISDOCID: <EP 0977054A2 I > 



EP0 977 054 A2 

the solution, take Z-transforms of the difference equation to get: 

4*0(z,u)=[(z + z + + u " 1 )]-4>(z,/y)+ P(z,u) ; (10) 

5 where: 

p i,H= m ijM - m h- i,a + mk ik - mk iM - 1 ; 

and: 

70 

<t>{z,u)=Z(0 jk ) (12) 

and 

rs P(z,u)=Z(p jk ) . 

Solving for <£<z.u) yields: 

20 *(z.u) = ¥^ - (13) 

20 4-[(z + z 1 ) + (" + u 1 )] 

Evaluating on the unit circle leads to: 



25 



DFT(pj H ) 

^* m,n~ 7 7 I3[ 7 TT 1 (14) 

30 and the solution is the first-quadrant portion of the inverse Discrete Fourier Transform (DFT), or: 

0 Jtk =(DFT(<t> mn y') Lk . (15) 

[0046] In this regard, it should be appreciated that the invention preferably uses mirror-image array dimensions which 

35 are powers of two so that the DFT can be evaluated using the Fast Fourier Transform (FFT). Moreover, a substantial 
computational advantage results from using many short FFT's rather than a single long one. As a result, decomposing 
the interferogram (e.g., 61) into small transform blocks (e.g., 64) is computationally more efficient than using large ones. 
[0047] Note that the denominator of Eq. (14) is zero when, and only when, m=n=0 . As a result, the solution is invalid 
if the value of the numerator is non-zero for m = n = 0 . Construction of the minor-image data, along with specified 

40 boundary conditions on the second differences, ensure that the second differences sum to zero over the transform 
block (e.g., 64). This results in a zero value for the (0,0) bin of the numerator transform, and hence a valid solution. 
[0048] It should also be noted that in a modified embodiment, step 70 in Fig. 1 may be completed after the phase 
unwrapping step 80. In such modified embodiment, the analysis of transform blocks (e.g., 64) to identify inconsistencies 
(step 70) may be completed in a manner illustrated in Fig. 3. Briefly, each transform block (e.g., 64) may be rewrapped 

45 to obtain rewrapped pixel data (step 71 ). The rewrapped pixel data may then be subtracted from the initial wrapped pixel 
data to obtain to obtain a difference function (step 72). The difference data may then be wrapped (step 73) and evalu- 
ated (step 74). Such evaluation may comprise a determination of whether the wrapped difference function is constant 
(step 74). If not, the corresponding data for an entire transform block (e.g.. 64) is discarded from further analysis in the 
process. If the difference function is constant, then the data corresponding with the evaluated transform block (e.g., 64) 

so is retained for further processing in the method. 

[0049] Returning now to the embodiment 1 0 of Fig. 1 , following the phase unwrapping step 80, integration constants 
for each of the retained transform blocks (e.g., 64) are determined (step 90). In this regard, it is noted that the phase 
values in the retained, unwrapped transform blocks (e.g., 64) are valid only to within an arbitrary constant. But since 
each transform block overlaps its neighbors (e.g., by a band of four samples), and since imaged terrain is unique, the 

55 average value of the unwrapped phase values in the overlap regions (e.g. . 65) must be the same for neighboring blocks, 
thereby allowing the constants for adjoining, retained blocks to be determined. 

[0050] More particularly, resolution of the constants may be achieved by first using a single known tie point to deter- 
mine the constant for one of the blocks. Then given the commonality of the overlap regions (e.g., 65), the constants for 
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alt of the other accessible, retained blocks (e.g., 64) may be determined by path-following. This procedure can be 
repeated for ail sets of disjointed transform blocks (i.e., sets that are not accessible via a path to a resolved transform 
block set due to the discarding of transform blocks in step 70 above), provided that a tie point is known with respect to 
each of the sets. The specific procedure is further described below with reference to Fig. 4. 
5 [0051] Let indices (u,v) count transform blocks (e.g., 64) comprising an interferogram (e.g., 61). Let L(u,v), R(u,v), 
U(u,v), and D(u,v) be the average unwrapped phase values in the left (decreasing v), right (increasing v), up (decreas- 
ing u), and down (increasing u) overlap regions (e.g., 65) of block (u.v). Let C(u,v) be the constant of integration lor block 
(u,v). Then, with respect to their four neighbors, the constants C(u,v) must satisfy: 



10 






Ru r v + Cu,v+ L u.v +1 + Cu,v +1 






15 






D u,v + Cu,v = Uu Cu+ 1,v 



Accordingly, if one knows a height value for a single data sample within a given transform block (u,v). i.e., a known "tie 
point", a corresponding value for a constant of integration for the transform block C(u,v) can be established (step 91). 

20 In turn, the above-noted expressions can be used to determine the constants for the neighboring blocks. More particu- 
larly, the unwrapped phase values can be averaged for the overwrap regions of the tie point transform block and for the 
corresponding overlap regions of adjacent transform blocks 92. Since the averages for corresponding overlap regions 
should be equal, the constant for the transform block containing the tie point can be used to determine the constant for 
adjacent transform blocks (step 93). In turn, a path-following approach can be used to determine the constants of inte- 

25 gration for all retained transform blocks (step 94). In the event that sets of retained transform blocks are inaccessible by 
path-following (e.g., as a result of the discarding of data for transform blocks in step 70), additional tie point information 
may be needed to determine the constants for the inaccessible transform blocks (see step 95). 
[0052] By way of example and referring to Fig. 2B, assume that the transform blocks 64 within region 66 are all dis- 
carded in step 70 due to inconsistencies found therewithin. In such a situation, region 67 of the interferogram portion 

30 62 may be rendered inaccessible via path-following. In such a situation, it would be necessary to have another known 
height to establish a tie point within the region 67 (step 96) in order to determine constants of integration for the trans- 
form blocks 64 with such region 67 (step 98). In the event that another tie point within region 67 is known, then path- 
following may be used within region 67 to determine integration constants for the various transform blocks with region 
67 (step 99), if not, the transform blocks within region 67 may be discarded from further use in the method. 

35 [0053] It is noted that, over the ensemble of blocks spanned by indices (u,v), the system of equations for the constants 
is overdetermined, making it theoretically possible for the resulting set of constants to be inconsistent which would 
result in erroneous height shears when the blocks are mosaiced together. Such inconsistencies do not occur in inter- 
ferograms from authentic terrain, but they can be induced if an interferogram (e.g., 61) is smoothed so severely that 
regions heavily tainted by noise are transformed into consistent phase data. This situation is minimized or avoided by 

40 ensuring that the span of the smoothing filter used in step 50 hereinabove is no more than about half of the size of a 
transform block. In this regard, the amount of smoothing in step 50 can be typically done with filter spans of 3-by-3, 5- 
by-5, or 7-by-7 data samples, all of which are compatible with the 17-by-17 or 33-by-33 transform blocks. 
[0054] After the constants of integration for retained transform blocks (e.g., 64) have been found, the unwrapped 
phase values may be converted to terrain height by using the constants and applying an interferometric phase-to-height 

45 conversion factor S (step 100), which is nominally constant over the span of transform block indices (u,v). Accordingly, 
the transform height values are given by: 

V aS -*/,*: (17) 

so where indices (j,k) span the entire interferogram (e.g., 61) with deletions for discarded blocks (i.e., resulting from step 
70). The value of S is a function of the specific IFSAR system and its collection geometry. Computing S is understood 
by those skilled in the art of radar interferometry. 

[0055] In many interferometric height-finding problems, multiple points with known terrain heights may be available 
and can be readily used in the process embodiment 10 to achieve enhanced accuracy (step 110). Since the height 
55 accuracies at these points may vary, they are incorporated by adjusting the determined heights within the retained 
transform blocks (e.g., 64) so as to minimize height differences in the overlap regions (e.g., 65). The mathematical for- 
mulation is a set of the difference equations which are iterated until convergence occurs. 

[0056] By way of explanation, let Ah(u,v) be the set of height adjustments required to incorporate the additional tie 
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points. Every pixel in transform block (u,v) is adjusted by Ah(u,v) for that block Each Ah(u,v) should satisfy the differ- 
ence equations: 



Ah uv + hL uv =Ah uv + hR u v . 1 



(18) 



Ah uv + hR uv =Ah uv +1 + hL uv _ n 



Ah u v + hU u v =Ah u v + hD u _ 1>v 



70 



Ah uv + hD uv =Ahj u +1 V +AA7 U W u+ 1iV , 



where hL, hR, hll, and hD are the average heights in the left, right, upper, and lower overlap regions of the transform 
blocks. Adding the equations and dividing by four yields: 



15 



(19) 



Ah =-■ 



I Ah + Ah , + Ah . + Ah , ) 



hL . , - hR 

u,v * 1 n.v 



hU 



V ♦ l.V 



20 



This equation is executed iteratrvely starting with each of the Ah(u,v) = 0 . During the execution, if a given transform 
block has been previously discarded (e.g., as containing inconsistent data in step 70), its indices are skipped. If a given 
block contains a known tie point, its Ah is pre-computed (i.e., to be equal to the difference between the known tie point 

25 value and the value determined by steps 20 through 1 00) and held constant during the iteration. The iteration continues 
until the sum of the squared Ah values changes negligibly from one iteration to the next (e.g., less than about 0.1%). 
after which. for each retained transform block, the corresponding Ah value is added to all pixels in the block. 
[0057] The Ah iterative solution described above serves to implement a weighted least-squares technique which min- 
imizes the sum of the squared height differences in the overlap regions (e.g., 65) of the retained blocks, and which is 

30 constrained by the known tie-point heights (i.e., to be equal to the difference between the known tie point value and the 
value determined by steps 20 through 100). 

[0058] A further effect of the Ah iteration is to smooth any inconsistencies which might arise in determining the con- 
stants of integration (step 90). For this reason, the iterative procedure may be executed even if only a single tie point is 
known. In this situation, only a single iteration is needed rf consistent constants of integration have been found. 
35 [0059] As will be appreciated, the foregoing description pertains to one illustrative embodiment of the present inven- 
tion. The various aspects of the described invention may be used in other embodiments and may be extended in other 
applications. All such extensions and applications are intended to be within the scope of the present invention. 



Claims 

40 

1. A method for determining height information from an interferogram of complex data samples generated from two 
complex images of an imaged terrain region acquired by an interferometric synthetic aperture radar system, com- 
prising: 



45 dividing said interferogram into a plurality of transform blocks, wherein each transform block comprises a plu- 

rality of complex data samples, and wherein at least a first of said transform blocks partially overlaps an adja- 
cent, second transform block in an overlapping region; 

unwrapping a phase difference value for each of the plurality of complex data samples comprising said first 
transform block in the overlapping region and for each comprising said second transform block in the overlap- 
so ping region; 

using a known height value corresponding with one of said plurality of complex data samples comprising said 
first transform block to establish an integration constant for the first transform block; and 
employing said integration constant for said first transform block and said unwrapped phase difference values 
to determine an integration constant for at least said second transform block. 

55 

2. A method as claimed in Claim 1 , said employing step comprising: 



averaging the unwrapped phase difference values for the overlapping region of the first transform block to 
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obtain a first average; 

averaging the unwrapped phase difference values for the overlapping region of the second transform block to 
obtain a second average; and 

using the first average and the second average to determine the integration constant for the second transform 
5 block. 

3. A method as claimed in Claim 1, further comprising: 

generating an integration constant for each of a first plurality of said transform blocks using said integration 
io constraints for said first and second transform blocks and path-following. 

4. A method as claimed in Claim 3, further comprising: 

smoothing the complex data samples comprising the interferogram using a spanning fitter of a predetermined 
75 size, wherein said predetermined size is less than about one-haH of said plurality of complex data samples 

comprising each said transform block. 

5. A method as claimed in Claim 1, wherein the total number of data samples comprising any one of said plurality of 
transform blocks is at least about 100. 

20 

6. A method as claimed in Claim 5, wherein each transform block comprises a matrix of o x p data samples, where o 
and p =» 2 n + 1 , and wherein n is an integer greater than 2. 

7. A method as claimed in Claim 5. wherein said overlapping region of comprises a total number of complex data 
25 samples that is between about 5% and 50% of the total number of data samples in the first transform block. 

8. A method as claimed in Claim 7, wherein said overlapping region is at least four complex data samples wide. 

9. A method as claimed in Claim 1, said unwrapping step comprising: 

30 

using a least-squares algorithm. 

10. A method as claimed in Claim 9. wherein said least squares algorithm is implemented by Fast Fourier Transform. 

35 1 1 . A method as claimed in Claim 1 , further comprising: 

analyzing each of said plurality of transform blocks separately to identify undesired data in the complex data 
samples comprising each transform block. 

40 12. A method as claimed in Claim 1 1 , wherein each of said complex data samples comprises a wrapped phase value, 
and wherein for each of said plurality of transform block(s) said analyzing step comprises: 

determining a difference value between the wrapped phase values of each set of adjacent complex data sam- 
ples comprising a transform block, 
45 comparing each said difference value to a predetermined value, wherein all complex data samples comprising 

a transform block are discarded from further use in the method when any single difference value corresponding 
with the transform block exceeds said predetermined value, and wherein all complex data samples comprising 
a transform block are retained for further use in the method when all difference values corresponding with the 
transform block are less than the predetermined value. 

50 

13. A method as claimed in Claim 12, wherein said predetermined value is less than about 135°. 

14. A method as claimed in Claim 12, wherein said predetermined value is less than about 90°. 

55 15. A method as claimed in Claim 12, further comprising: 

generating an integration constant for each of a first plurality of said retained transform blocks using said inte- 
gration constants for said first and second transform blocks and path-following; 
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using a known height value corresponding with one of said plurality of complex data samples comprising a third 
transform block to establish an integration constant for the third transform block, wherein said third transform 
block is separated from said first plurality of said retained transform blocks; 

generating an integration constant for each of a second plurality of said retained transform blocks using said 
s integration constant for said third transform block and path-following. 

16. A method as claimed in Claim 1 , further comprising: 

using said integration constant for said second transform block to determine height values for the data samples 
10 comprising the second transform block; and 

employing a known height value corresponding with one of said plurality of complex data samples comprising 
one of said transform blocks, other than said first transform block, to determine a height adjustment value for 
application to said heights determined for said complex data samples comprising said second transform block. 

15 
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(54) Determining height information 

(57) A method is provided for phase unwrapping in 
coherent imaging systems which use two complex 
images to form an interferogram. In the method, the 
interferogram of phase differences is divided into a rela- 
tively large plurality of transform blocks, wherein each 
transform block comprises a matrix of o x p data sam- 
ples, with o and p = 2 n + 1, and n being an integer 
greater than 2. Further, adjacent transform blocks are 
preferably defined to partially overlap. The wrapped 
phase values of adjacent data samples within each 
transform block are then compared and if the difference 
therebetween exceeds a predetermined value, the 
entire block of data samples may be discarded for fur- 
ther use. Data samples of the retained transform blocks 
are then unwrapped to obtain phase values via an 
unweighted, least-squares technique implemented by 
Fast Fourier Transform. Using a known height value cor- 
responding with a single data sample, an integration 
constant can be determined for a corresponding first 
transform block. Path -following may then be used to 
determine an integration constant for the other retained 
transform blocks. The phase values and integration con- 
stants may then be employed to determine height val- 
ues for the interferogram. The method yields enhanced 
accuracy with reduced computational burden. 
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